##### for Figure 4 b,d

rm(list=ls())
options(stringsAsFactors = F)
library(Seurat)
library(ggplot2)
library(dplyr)
library(clustree)
#setwd("D:/OneDrive - zju.edu.cn/data/public_database/scRNA/E-MTAB-7703")
setwd("/Users/gin/Documents/onedrive/data/public_database/scRNA/E-MTAB-7703/")
##### Data format ##### ---------
sce <- CreateSeuratObject(Read10X('./sce/'),"sce")


##### reduction ##### ---------
library(Seurat)

sce <- NormalizeData(sce) %>% 
  FindVariableFeatures(selection.method = "vst",nfeatures = 3000) %>% 
  ScaleData() %>% 
  RunPCA(npcs = 30, verbose = T)

sce <- FindNeighbors(sce, dims = 1:30)
sce <- FindClusters(object = sce,resolution = c(seq(0,1,.1)))
clustree(sce)
sce <- FindClusters(sce, resolution = 0.2)
# umap/tsne降维
sce <- RunTSNE(sce, dims = 1:20)
sce <- RunUMAP(sce, dims = 1:20)
DimPlot(sce,reduction = "tsne",pt.size = 1,label=T)
DimPlot(sce,reduction = "umap")
FeaturePlot(sce,"Grem1",reduction = "tsne",order=T,pt.size = 1)
FeaturePlot(sce,"Ccl19",reduction = "tsne",order=T,pt.size = 1)

##### stackvlnplot ##### -----
library(Seurat)
library(ggplot2)
### Grem1+ and Grem1- ###
g_pos = subset(sce,ident = 3)
g_neg = subset(sce,ident = c(0,1,2,4:11))

meta_gpos = as.data.frame(colnames(g_pos))
colnames(meta_gpos) = "cell"
meta_gpos$group="g_pos"

meta_gneg = as.data.frame(colnames(g_neg))
colnames(meta_gneg) = "cell"
meta_gneg$group="g_neg"

meta = rbind(meta_gpos,meta_gneg)
meta_sce = as.data.frame(colnames(sce))
colnames(meta_sce) = "cell"

meta = merge(meta_sce,meta,by="cell",sort = F)
rownames(meta) = meta[,1]
#meta = meta[,-1]

sce <- AddMetaData(object = sce, metadata = meta[,2], col.name = "group")


### 运行
#genelist = c('Grem1',"Ccl21a",
#             "Tcf21","Wt1","Tlx1","Ly6a",
#             "Dpt","Gsn","Thy1","Il33",
#             "Madcam1","Cxcl13","Igfbp3")
genelist = c('Grem1',
             "Tcf21","Tlx1",
             "Bst1","Ccl19","Ccl21a","Crym","Il33",
             "Madcam1","Cxcl13","Igfbp3")
VlnPlot(sce, features = genelist, pt.size = 0 ,stack=T,group.by = "group") 

#####


rm(list=ls())
options(stringsAsFactors = F)
library(Seurat)
library(ggplot2)
#setwd("D:/OneDrive - zju.edu.cn/data/public_database/scRNA/E-MTAB-7703")
setwd("/Users/gin/Documents/onedrive/data/public_database/scRNA/E-MTAB-7703/")
##### Data format ##### ---------
sce <- CreateSeuratObject(Read10X('./sce/'),"sce")


##### reduction ##### ---------
library(Seurat)

sce <- NormalizeData(sce) %>% 
  FindVariableFeatures(selection.method = "vst",nfeatures = 3000) %>% 
  ScaleData() %>% 
  RunPCA(npcs = 30, verbose = T)

sce <- FindNeighbors(sce, dims = 1:30)
sce <- FindClusters(object = sce,resolution = c(seq(0,1,.1)))
clustree(sce)
sce <- FindClusters(sce, resolution = 0.2)
# umap/tsne降维
sce <- RunTSNE(sce, dims = 1:20)
sce <- RunUMAP(sce, dims = 1:20)
DimPlot(sce,reduction = "tsne",pt.size = 1,label=T)
DimPlot(sce,reduction = "umap")
FeaturePlot(sce,"Grem1",reduction = "tsne",order=T,pt.size = 1)
FeaturePlot(sce,"Ccl19",reduction = "tsne",order=T,pt.size = 1)

##### stackvlnplot ##### -----
library(Seurat)
library(ggplot2)
### Grem1+ and Grem1- ###
g_pos = subset(sce,ident = 3)
g_neg = subset(sce,ident = c(0,1,2,4:11))

meta_gpos = as.data.frame(colnames(g_pos))
colnames(meta_gpos) = "cell"
meta_gpos$group="g_pos"

meta_gneg = as.data.frame(colnames(g_neg))
colnames(meta_gneg) = "cell"
meta_gneg$group="g_neg"

meta = rbind(meta_gpos,meta_gneg)
meta_sce = as.data.frame(colnames(sce))
colnames(meta_sce) = "cell"

meta = merge(meta_sce,meta,by="cell",sort = F)
rownames(meta) = meta[,1]
#meta = meta[,-1]

sce <- AddMetaData(object = sce, metadata = meta[,2], col.name = "group")


### 运行
#genelist = c('Grem1',"Ccl21a",
#             "Tcf21","Wt1","Tlx1","Ly6a",
#             "Dpt","Gsn","Thy1","Il33",
#             "Madcam1","Cxcl13","Igfbp3")
genelist = c('Grem1',
             "Tcf21","Tlx1",
             "Bst1","Ccl19","Ccl21a","Crym","Il33",
             "Madcam1","Cxcl13","Igfbp3")
VlnPlot(sce, features = genelist, pt.size = 0 ,stack=T,group.by = "group") 

#####